Efficient  greedy  algorithms  for  high-dimensional  parameter 
spaces  with  applications  to  empirical  interpolation  and  reduced 

basis  methods 

Jan  S.  Hesthaven  *  Benjamin  Stamm  '  Shun  Zhang  * 

September  9,  2011 

Abstract.  We  propose  two  new  and  enhanced  algorithms  for  greedy  sampling  of  high¬ 
dimensional  functions.  While  the  techniques  have  a  substantial  degree  of  generality,  we 
frame  the  discussion  in  the  context  of  methods  for  empirical  interpolation  and  the  devel¬ 
opment  of  reduced  basis  techniques  for  high-dimensional  parametrized  functions.  The  first 
algorithm,  based  on  a  assumption  of  saturation  of  error  in  the  greedy  algorithm,  is  shown 
to  result  in  a  significant  reduction  of  the  workload  over  the  standard  greedy  algorithm.  In 
an  improved  approach,  this  is  combined  with  an  algorithm  in  which  the  train  set  for  the 
greedy  approach  is  adaptively  sparsefied  and  enriched.  A  safety  check  step  is  added  at 
the  end  of  the  algorithm  to  certify  the  quality  of  the  basis  set.  Both  these  techniques  are 
applicable  to  high-dimensional  problems  and  we  shall  demonstrate  their  performance  on  a 
number  of  numerical  examples. 


1  Introduction 

Approximation  of  a  function  is  a  generic  problem  in  mathematical  and  numerical  analysis, 
involving  the  choice  of  some  suitable  representation  of  the  function  and  a  statement  about  how 
this  representation  should  approximate  the  function.  A  traditional  approach  is  polynomial 
representation  where  the  polynomials  coefficients  are  chosen  to  ensure  that  the  approximation  is 
exact  at  certain  specific  points,  recognized  as  the  Lagrange  form  of  the  interpolating  polynomial 
representation.  Such  representations,  known  as  linear  approximations,  are  independent  of  the 
function  being  approximated  and  have  been  used  widely  for  centuries.  However,  as  problems 
becomes  complex  and  high-dinrensional,  the  direct  extension  of  such  ideas  quickly  becomes 
prohibitive. 

More  recently,  there  has  been  an  increasing  interest  in  the  development  of  nonlinear  ap¬ 
proximations  in  which  case  the  approximation  is  constructed  in  a  problem  specific  manner  to 
reduce  the  overall  computational  complexity  of  constructing  and  evaluating  the  approximation 
to  a  given  accuracy.  In  this  setting,  the  key  question  becomes  how  to  add  an  element  to  the 
existing  approximation  such  that  the  new  enriched  approximation  improves  as  much  as  possi¬ 
ble,  measured  in  some  reasonable  manner.  This  approach,  known  as  a  greedy  approximation, 
seeks  to  maximize  a  given  function,  say  the  maximum  error,  and  enrich  the  basis  to  eliminate 
this  specific  error,  hence  increasing  the  accuracy  in  an  optimal  manner  when  measured  in  the 
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maximum  norm.  Such  a  greedy  approach  has  proven  themselves  to  be  particularly  valuable  for 
the  approximation  of  high-dinrensional  problems  where  simple  approaches  are  excluded  due  to 
the  curse  of  dimension.  For  a  detailed  recent  overview  of  such  ideas  in  a  general  context,  we 
refer  to  [15] . 

In  this  work  we  consider  greedy  algorithms  and  improvements  of  particular  relevance  to 
high-dinrensional  problems.  While  the  ideas  are  of  a  general  nature,  we  motivate  and  frame  the 
discussion  in  the  context  of  reduced  basis  methods  (RBM)  and  empirical  interpolation  methods 
(EIM)  in  which  the  greedy  approximation  approach  plays  a  key  role.  In  the  generic  greedy 
approach,  one  typically  needs  a  fine  enough  train  set  Strain  C  T>  over  which  a  functional  has  to 
be  evaluated  to  select  the  next  element  of  the  approximation.  When  the  number  of  parameters 
is  high,  the  size  of  this  train  set  quickly  becomes  considerable,  rendering  the  computational 
cost  substantial  and  perhaps  even  prohibitive.  As  consequence,  since  a  fine  enough  train  set 
is  not  realistic  in  practice,  one  is  faced  with  the  problem  of  ensuring  the  quality  of  the  basis 
set  under  a  non-rich  enough  train  set.  It  is  worth  noting  that  when  dealing  with  certain  high 
dimensional  problems,  one  may  encounter  the  situation  that  the  optimal  basis  set  itself  is  of 
large  size.  This  situation  is,  however,  caused  by  the  general  complexity  of  the  problems  and 
we  shall  not  discuss  this  further.  Strategies  for  such  cases  are  discussed  in  see  [6,  7]. 

In  this  paper,  we  propose  two  enhanced  greedy  algorithms  related  to  the  search/loop  over 
the  large  train  set.  The  first  algorithm  utilizes  a  saturation  assumption,  based  on  the  assump¬ 
tion  that  the  greedy  algorithm  converges,  i.e.,  with  enough  bases,  the  error  will  decrease  to 
zero.  It  is  then  reasonable  to  assume  that  the  error  (or  the  error  estimator  in  the  case  of  the 
reduced  basis  method)  is  likewise  decreasing  in  some  sense.  With  a  simple  and  reasonable 
saturation  assumption  on  the  error  or  the  error  estimator,  we  demonstrate  how  to  modify  the 
greedy  algorithm  such  only  errors/error  estimators  are  computed  for  those  points  in  Strain 
with  a  large  enough  predicted  error.  With  this  simple  modification,  the  total  workload  of  the 
standard  greedy  algorithm  is  significantly  reduced. 

The  second  algorithm  is  an  adaptively  enriching  greedy  algorithm.  In  this  approach,  the 
samples  in  the  train  set  is  adaptively  removed  and  enriched,  and  a  safety  check  step  is  added 
at  the  end  of  the  algorithm  to  ensure  the  quality  of  the  basis  set.  On  each  step  of  searching 
a  new  parameter  for  a  basis,  the  size  of  the  train  set  is  maintained  at  a  reasonable  number. 
This  algorithm  can  be  applied  to  problems  with  high  number  of  parameters  with  substantial 
savings. 

What  remains  of  this  paper  is  organized  as  follows.  In  Section  2,  we  discuss  the  role  greedy 
sampling  plans  in  different  computational  methods,  exemplified  by  empirical  interpolation  and 
reduced  basis  methods,  to  highlight  shortcomings  of  a  naive  approach  and  motivate  the  need 
for  improved  methods.  This  sets  the  stage  for  Section  3  where  we  discuss  the  details  of  two 
enhanced  greedy  techniques.  This  is  followed  in  Section  4  and  5  by  a  number  of  detailed 
numerical  examples  for  the  empirical  interpolation  methods  and  reduced  basis  techniques, 
respectively,  to  illustrate  the  advantages  of  using  these  new  methods  for  problems  with  both 
low  and  high-dinrensional  parameter  spaces.  Section  6  contains  a  few  concluding  remarks. 

2  On  the  need  for  improved  greedy  methods 

In  the  following  we  give  a  brief  background  on  two  different  computational  techniques,  both 
of  which  rely  on  greedy  approximation  techniques,  to  serve  as  motivation  for  the  subsequent 
discussion  of  the  new  greedy  techniques. 
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2.1  Reduced  basis  methods 


Many  applications  related  to  computational  optimization,  control,  and  design  require  the  abil¬ 
ity  to  rapidly  and  accurately  solve  parameterized  problems  many  times  for  different  parameter 
values  within  a  given  parametric  domain  PcP.  While  there  are  many  suitable  methods  for 
this,  we  shall  focus  here  on  the  reduced  basis  method  (RBM)  [12,  14]  which  has  proven  itself 
to  be  an  very  accurate  and  efficient  method  for  such  scenarios. 

For  any  /x  £  T>,  the  goal  is  to  evaluate  an  output  functional  s(/x)  =  £(tt(/x); /x),  where 
xt(/x)  el  is  the  solution  of 

a(u(n),v,  /x)  =  f(v,  /x),  Vwel  (2.1) 


for  some  parameter  dependent  bilinear  and  linear  forms  a  and  /  and  X  is  a  suitable  functional 
space. 

Let  X*e  be  a  finite  element  discretization  subspace  of  X.  Here,  finite  elements  are  used  for 
simplicity,  and  other  types  of  discretizations  can  likewise  be  considered.  For  a  fixed  parameter 
/x  E  T>,  let  E  Xfe  be  the  numerical  solution  of  the  following  Galerkin  problem, 

a(w/e(/x),  v,  /x)  =  f(v;  /*),  Vu  G  Xfe,  (2.2) 

and  let  s^e{n)  =  £(xHe(/x);  /x)  be  the  corresponding  output  functional  of  interest. 

Both  the  variational  problem  (2.1)  and  the  approximation  problem  (2.2)  are  assumed  to 
be  well-posed.  The  following  inf-sup  stabilities  are  satisfied  for  /x-dependent  positive  constants 
/3(/x)  and  respectively: 

/?(/x)  =  inf  sup  and  /3/e(/x)  =  inf  sup  ,  (2.3) 

U^X  V£X  IMU'IMIx  U&Xfe  v&Xfe  IMIX-fe  IMIWe 


where  ||  •  |  x  and  ||  •  || Xfe  are  norms  of  the  spaces  X  and  X*e,  respectively. 

For  a  collection  of  N  parameters  Sx  =  {/x1,  •  •  •  ,  /x^}  in  the  parameter  domain  T>  C  Mp, 
let  Wx  =  {u-^/x1),  •  •  •  ,  u^e(fiN)},  where  u^e(/x*)  is  the  numerical  solution  of  problem  (2.2) 
corresponding  to  the  parameter  values  /x*,  for  1  <  i  <  N.  Then,  define  the  reduced  basis  space 
as  X$  =  span{T'Rv}- 

The  reduced  basis  approximation  is  now  defined  as:  For  a  /x  G  V,  find  tt)y(/x)  G  Xx  such 
that 

«(«#(/*)>*>;  aO  =  f(v]  /*),  Vv  £  xtf,  (2.4) 

with  the  corresponding  value  of  the  output  functional 

$(/*)  =  (2.5) 


Define  the  error  function  e(/x)  =  urx(fj,)  —  tFe(/x)  G  X*e  as  the  difference  between  the 
reduced  basis  (RB)  solution  ux(fi)  and  the  highly  accurate  finite  element  solution  iHe(/x). 
The  residual  r(v;  /x)  G  (X^e)'  is  defined  as 

r(v;  /x)  :=  f{v;  /x)  -  a(4,r;g),  Vu  G  Xfe,  (2.6) 


and  its  norm  as 


II^saOII (xf*y  ■=  SUP 

V&Xfe 

We  define  the  relative  estimator  for  the  output  as 


r(v;p) 

IMLy /«’ 


?y(/x,  WN) 


V)\\(xfey¥fe(-,  riWjXfey 
/5/e(/x)|s^(/x)| 


(2.7) 


(2.8) 
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Other  types  of  error  estimators  can  also  be  used,  see  e.g.,  [14]. 

To  build  the  parameter  set  Sn,  the  corresponding  basis  set  Wn  and  the  reduced  basis  space 
X$,  a  greedy  algorithm  is  used.  For  a  train  set  retrain  C  T>,  which  consists  of  a  fine  discretiza¬ 
tion  of  V  of  finite  cardinality,  we  first  pick  a  /i1  G  Strain,  and  compute  the  corresponding  basis 
ufe(nx).  Let  S\  =  {/x1},  W\  =  and  X[b  =  span {r/e(/x1)}.  Now,  suppose  that  we 

already  found  N  points  in  Strain  to  form  Sn,  the  corresponding  Wn  and  X$,  for  some  integer 
N  >  1.  Then,  choose 

l^N+1  ■=  argmaxMeStro,n??(M;  WN),  (2.9) 

to  fix  the  next  sample  point  and  let  Sn+i  '■=  Sn  U  {/^.A  +1 } .  We  then  build  the  corresponding 
spaces  Wn+i  and  Xy^‘{.  The  above  procedure  is  repeated  until  N  is  large  enough  that 
max/xg=train  r]{n\  Wn)  is  less  than  a  prescribed  tolerance. 

For  this  approach  to  be  successful,  it  is  essential  that  the  training  set  is  sufficiently  fine, 
i.e.,  for  problems  with  many  parameters,  the  size  of  the  train  set  Strain  could  become  very 
large.  Even  with  a  rapid  approach  for  evaluating  r/(/i;  Wn)  for  all  fi  G  Strain  the  cost  of  this 
quickly  becomes  a  bottleneck  in  the  construction  of  the  reduced  basis. 

2.2  Empirical  interpolation  method 

One  of  the  main  attractions  of  the  reduced  basis  method  discussed  above  becomes  apparent  if 
we  assume  that  the  parameter  dependent  problem  (2.1)  satisfies  an  affine  assumption,  that  is, 


Q  a  Qf  Qi 

a(u,  v-  p)  =  0“(/x)a*(u,  v),  f(v,  /x)  =  &{ (n)fi(v),  and  £(v,  n)  =  ^  © f(pK;(u), 

i=  1  i=  1  i= 1 

(2.10) 

where  0“,  (~)j,  and  @(  are  /x-dependent  functions,  and  cq,  /) ,  are  p-in dependent  forms. 
With  this  assumption,  for  a  reduced  basis  space  X$  with  N  bases,  we  can  applu  the  so-called 
offline/online  strategy.  In  the  offline  step,  one  can  precompute  the  matrices  and  vectors  related 
to  forms  a,;,  /*,  and  li,  for  i  =  1,  •  •  •  ,  N.  The  cost  of  this  may  be  substantial  but  is  done  only 
once.  In  the  online  step,  we  now  construct  the  matrices  and  vectors  in  the  reduced  basis 
formulation  (2.4),  solve  the  resulting  reduced  basis  problem,  and  then  evaluate  the  output 
functional  (2.5).  The  amount  of  work  of  the  online  step  is  independent  of  the  degrees  of 
freedom  of  X^e,  and  only  depends  on  the  size  of  reduced  basis  N  and  the  affine  constants 
Qa ,  Qf,  and  Qg.  Hence,  for  a  fixed  p  G  V,  the  computation  77(7x5  Wn)  includes  the  solving 
procedure  of  the  reduced  basis  problem,  the  evaluation  of  the  residual  (and  output  functional) 
in  the  dual  norm,  and  a  possible  linear  programming  problem  to  evaluate  see  [2].  As 

already  anticipated,  the  amount  of  work  does  not  depend  on  the  size  of  X*e,  but  only  on  N 
and  is,  hence,  very  fast. 

However,  when  the  parameter  dependent  problem  does  not  satisfy  the  affine  assumption 
(2.10),  this  key  benefit  is  lost.  To  circumvent  this,  the  empirical  interpolation  method  (EIM) 
[2,  9,  8]  has  been  developed  to  enable  one  to  treat  the  non- affine  operators  and  approximate 
them  on  the  form  (2.10)  to  maintain  computational  efficiency. 

To  explain  the  EIM,  consider  a  parameter  dependent  function  T  :  H  x  D  -+  1  or  T  : 
H  x  V  — >  C.  The  EIM  is  introduced  in  [2,  9,  11]  and  serves  to  provide  parameter  values 
Sn  =  {h1)  •  -  ■  5  t^N}  such  that  the  interpolant 


N 

XatC^Xx;^)  :=  ^  (3n(fj,)qn(x) 

71=  1 


(2.11) 
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is  an  accurate  approximation  to  ^(x;  fi)  on  H  x  V. 

The  sample  points  Sn  are  chosen  by  the  following  greedy  algorithm.  Again,  using  a 
train  set  Strain  C  T>  which  consists  of  a  fine  discretization  of  V  of  finite  cardinality,  we 
first  pick  a  /x1  6  Strain-,  compute  x1  =  argnraxxef2  B(x;  /x1)  and  the  corresponding  basis 
<71  (-)  =  /i1)/iF(x1;  /i1).  Then,  let  S\  =  {/x1}  and  1-Ti  =  [q\ }. 

Now,  suppose  that  we  already  found  N  points  in  Strain  to  form  Sn  and  Wn  =  {q\ , . . . ,  qw} 
such  that  span{(/i, . . . ,  qjy}  =  span{Jr(-;  /x1), . . . ,  for  some  integer  N  >  1.  We  further 

assume  that  a  set  of  N  points  T/v  =  {x1, . . . ,  xA  }  is  given  such  that 

1  if  *  =  j, 

0  if  i  <  j ,  (2.12) 

q-j(xl)  otherwise, 

for  all  i ,  j  =  1, . . .  ,N.  Then,  we  denote  the  lower  triangular  interpolation  matrix  B ?A  =  c/j(x*), 
z,  j  =  1, . . . ,  N,  to  define  the  coefficients  {f3n(^)}^=1,  for  a  given  fx  6  P,  which  are  the  solution 
of  the  linear  system 

N 

=  ^r(xi;^),  V*  =  i,  . . . ,  iv. 

i=i 

The  approximation  of  level  N  of  (x)  is  given  by  the  interpolant  defined  by  (2.11).  We  then 
set 

r](/j,;WN)  ■=  \\B(-;ix)  -  /x)  ||L°°(f2) 

and  choose 

VN+1  ■=  argmaxMgStra,n??(/x;  WN ),  (2.13) 

to  fix  the  next  sample  point  and  let  SW+i  :=  Sn  U  { /xA  “t_1 } .  The  interpolation  point  xjV+1  is 
defined  by 

x^1  =  argmaxxgn(jF(.;/x)  -  2jv(^)(-;  /*)) 
and  the  next  basis  function  by 

QN+1  ^(x^+i;  /x)  -  Tjv(^7) (xAr+1;  /x) ' 

By  construction,  we  therefore  satisfy  the  condition  (2.12)  since  the  interpolation  is  exact  for  all 
previous  sample  points  in  Sn-  The  algorithm  is  stopped  once  the  error  maxM6strai„  r/(/x;  Wv) 
is  below  some  prescribed  tolerance.  As  one  can  observe,  the  EIM  also  uses  a  greedy  algo¬ 
rithm  to  choose  the  sample  points  with  only  slight  differences  to  the  case  of  the  reduced  basis 
method.  Hence,  in  the  case  of  a  high  dimensional  parameter  space,  the  computational  cost  of 
constructing  the  empirical  interpolation  can  be  substantial. 

3  Improved  greedy  algorithms 

Having  realized  the  key  role  the  greedy  approach  plays  in  the  two  methods  discussed  above, 
it  is  clear  that  even  if  the  greedy  approach  is  used  only  in  the  offline  phase,  it  can  result  in  a 
very  considerable  computational  cost,  in  particular  in  the  case  of  a  high-dinrensional  parameter 
space.  Let  us  in  the  following  discuss  two  ideas  aimed  to  help  reduce  the  computational  of  the 
greedy  approach  in  this  case. 


ffi(x'0  =  S 
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3.1  A  typical  greedy  algorithm 

To  make  the  algorithm  more  general  than  discussed  in  the  above,  we  make  several  assumptions. 
For  each  parameter  /x  in  the  parameter  domain  V  C  Mp,  a  /x-dependent  basis  function  u(/x) 
can  be  computed.  Let 

SN  =  {/*V-  -  ,HN} 

be  a  collection  of  N  parameters  in  T>  and 

WN  =  M/x1),---  ,v(nN)} 

be  the  collection  of  N  basis  functions  based  on  the  parameter  set  S n-  For  each  parameter 
/x  £  V,  suppose  that  we  can  compute  an  error  estimator  r/(/x;  W n)  of  the  approximation  based 
on  Wn- 

The  following  represent  a  typical  greedy  algorithm. 


Input:  A  train  set  Strain  C  V,  a  tolerance  tol  >  0 
Output:  Sn  and  Wn 

1:  Initialization:  Choose  an  initial  parameter  value  /x1  £  Strain,  set  Si  =  {/x1},  computj 
^(/x1),  set  W\  =  {^(/x1)},  and  N  =  1  ; 
while  max^gs  r/(/x;  Wn)  >  tol  do 
For  all  /x  £  Strain,  compute  r/(/x;  Wat)  ; 

Choose  /x^1  =  argmaxMeS(ro.nr/(^;  Wat); 

Set  5,v+i  =  Sjv  U  {/x^1}; 

Compute  v(/xA'+1),  and  set  VFat+i  =  Wn  U  (u(/xiV+1)}; 

Ni-N  +  1] 

end  while 


Algorithm  1:  A  Typical  Greedy  Algorithm 


Note  that  Sn  and  Wn  are  hierarchical: 

SN  CSM,  WNCWM  if  1  <  N  <  M. 


3.2  An  improved  greedy  algorithm  based  on  a  saturation  assumption 

As  mentioned  before,  on  step  3  of  the  greedy  algorithm  1,  we  have  to  compute  ??(/x;  Wn)  for 
every  /x  £  Strain ■  When  the  size  of  Strain  is  large,  the  computational  cost  of  this  task  is  very 
high.  Fortunately,  in  many  cases,  the  following  saturation  assumption  holds: 

Definition  3.1.  Saturation  Assumption 

Assume  that  r/(/x;  Wn)  >  0  is  an  error  estimator  depending  on  a  parameter  /x  and  a  hier¬ 
archical  basis  space  Wn-  If  the  following  property 

??(/x;  Wm)  <  Csa  r/(/x;  Wn)  for  some  Csa  >  0  for  all  0  <  N  <  M  (3.14) 

holds,  we  say  that  the  Saturation  Assumption  is  satisfied. 

Remark  3.2.  When  Csa  =  1,  it  implies  that  rj{pL\  Wn)  is  not  increasing  for  a  fixed  /x  and 
increasing  N.  When  Csa  <  1,  this  is  a  more  aggressive  assumption,  ensuring  that  t](ij,;Wn) 
is  strictly  decreasing.  This  assumption  with  Csa  <  1  is  very  common  in  the  adaptive  finite 
element  method  community,  see  [1],  The  assumption  Csa  >  1  is  a  more  relaxed  assumption, 
allowing  that  rj^pi ;  Wn)  might  not  be  monotonically  decreasing,  but  can  oscillating.  Since  the 
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underlying  assumption  of  the  greedy  algorithm  is  that  g(pi,  Wn)  will  converge  to  zero  as  N 
approaches  infinity,  we  can  safely  assume  that  even  if  ?7(/x,  Wn)  might  exhibit  intermittend 
non-monotone  as  N  is  increasing,  overall  it  is  decreasing. 

Utilizing  this  assumption,  we  can  design  an  improved  greedy  algorithm.  First,  for  each 
parameter  value  pi  E  Strain,  we  create  an  error  profile  r]saVed( /■*)■  For  instance,  initially,  we  can 
set  ijSaved{ t1)  =  hil1',  Wo)  =  oo.  Now  suppose  that  Sn  and  Wn  are  determined  and  that  we 
aim  to  find  the  next  sample  point  piN+1  =  argmax/xgStra.j r/(/i;  Wn).  When  pi  runs  through 
over  the  train  set  Strain ,  we  naturally  keep  updating  a  temporary  maximum  (over  Strain), 
until  we  have  searched  the  whole  set  Strain-  In  this  loop,  we  may  require  that,  for  each  pi, 
Vsaved(^)  =  viw,  W l)  for  some  L  <  N.  Now,  since  the  Saturation  Assumption  77(7x5  Wn)  < 
Csa  77(7x5  Wl)  for  L  <  N  holds  and  if  Csa  rjsaved(^)  is  less  than  the  current  temporary  maximum, 
77(7X5  Wn)  can  not  be  greater  than  the  current  temporary  maximum.  Hence,  we  may  skip  the 
computation  of  77(7X5  Wn),  and  leave  risaved(pi)  untouched.  On  the  other  hand,  if  Csagsavedin) 
is  greater  than  the  current  temporary  maximum,  it  is  potentially  the  maximizer.  Hence,  we 
compute  77(/x,  Wn),  update  rjsavedil1),  and  compare  it  with  the  current  maximum  to  see  if  an 
update  of  the  current  temporary  maximum  is  needed.  We  notice  that  if  we  proceed  in  this 
manner,  then  the  requirement  that  for  each  pi,  r]saved{hL)  =  17(7x5  Wl)  for  some  L  <  N  holds. 

The  algorithm  2  in  pseudo-code  of  the  Saturation  Assumption  based  algorithm  is  given  as: 


Input:  A  train  set  Strain  C  V,  Csa,  and  a  tolerance  tol 
Output:  Sn  and  Wn 

1:  Choose  an  initial  parameter  value  pi1  E  Strain,  set  S\  =  {pi1}',  compute  v(pi1),  set 
Wi  =  {^(/x1)},  and  IV  =  1; 

2:  Set  a  vector  r)saVed  with  r)saVed(lA  =  00  for  all  pi  ^  Strain  5 
3:  while  maxMeStrojn  gsaVed{^)  >  tol  do 
4:  errOVtmpmax  —  0) 

5:  for  all  pi  E  Strain  do 

6:  if  Csar]saved{n)  errortmpmax  then 

7:  Compute  77(7x5  WN)  ,  and  let  T]saved{pi)  =  r/(pi,  WN)', 

8:  if  risaved(v)  >  errortmpmax  then 

9:  errortmpmax  —  hsavedil-^) ,  and  let  pimax  —  U- 

10:  end  if 

11:  end  if 

12:  end  for 

13:  Choose  piN+1  =  pimax,  set  SN+1  =  Sn  u  {piN+1}', 

14:  Compute  v(piN+1),  set  Wn+ 1  =  Wn  U  {v(piN+1)}; 

15:  N  <—  N  +  1; 

16:  end  while 

Algorithm  2:  A  greedy  algorithm  based  on  a  saturation  assumption 


Remark  3.3.  Initially,  we  set  rjsavedil1)  =  00  for  any  pi  E  Strain  to  make  ensure  every  ifipi,  Wfi) 
will  be  computed. 

Remark  3.4.  Due  to  round-off  errors,  the  error  estimator  may  stagnate  even  if  we  add  more 
bases,  or  the  greedy  algorithm  will  select  some  point  already  in  Sn-  In  this  case,  the  greedy 
algorithm  should  be  terminated. 
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3.3  An  adaptively  enriching  greedy  algorithm 

Even  though  the  above  saturation  assumption  based  algorithm  has  the  potential  to  reduce 
the  overall  computational  cost,  it  may  still  be  too  costly  for  problems  with  a  high  number  of 
parameters.  Notice  that,  in  the  above  algorithms,  the  assumption  that  the  train  set  Strain  is  a 
fine  enough  discrete  subset  of  V  is  essential;  otherwise,  we  might  miss  some  phenomena  that 
are  not  represented  by  E '.train-  The  consequence  of  this  is  that  for  the  final  sets  of  Sn  and  Wn, 
there  may  be  some  parameter  p,  mV  but  not  in  retrain  such  that  ififi,  Wn)  is  greater  (or  even 
much  greater)  than  tol. 

Thus,  for  high  dimensional  parameter  problems,  we  will  likely  have  to  face  the  problem 
that  the  train  set  is  not  rich  enough.  To  address  this  problem,  we  first  build  the  set  of  basis 
Wn  based  on  a  train  set  with  a  relatively  small  number  of  points. 

To  ensure  that  the  set  of  basis  functions  corresponding  to  this  train  set  is  good  enough, 
we  add  a  ” safety  check”  step  at  the  end,  that  is,  we  test  the  bases  by  a  larger  number  of  test 
parameters,  to  see  if  the  resulting  error  estimators  are  less  than  the  prescribed  tolerance  on  this 
larger  parameter  set  too.  If  for  all  test  points,  the  estimated  errors  are  less  then  the  tolerance, 
we  may  conclude  that  the  basis  set  is  good  enough.  Otherwise,  we  add  the  first  failed  test 
point  (whose  error  is  larger  than  the  tolerance)  to  Sn,  and  redo  the  ” safety  check”  step  until 
the  set  of  basis  Wn  passes  the  ” safety  check”. 

For  problems  with  a  high  number  of  of  parameters,  it  is  in  practice  hard  to  construct  a  rich 
enough  train  set.  First,  it  is  almost  impossible  to  construct  a  tensor  product  based  train  set 
for  p  >  10.  Even  if  we  only  put  3  points  for  each  parameter,  which  is  of  course  far  from  rich, 
a  train  set  of  311  =  177’147  points  is  quite  big.  For  a  bigger  p,  the  curse  of  dimensionality  is 
inevitable. 

Instead  of  starting  from  a  tensor  product  train  set,  we  consider  a  (quasi-) random  train 
set.  However,  the  random  train  set  faces  the  same  problem:  it  is  far  from  ’’rich  enough”  for 
a  high  dimensional  parameter  problem.  Fortunately,  we  can  adaptively  change  the  train  set 
by  removing  useless  points  and  enriching  new  points.  Notice  that,  after  we  have  determined 
a  set  of  basis  functions,  the  estimated  errors  corresponding  to  some  points  in  the  train  set  are 
already  smaller  than  the  prescribed  tolerance.  There  is  no  value  in  retaining  these  points  in 
Sn  they  should  be  removed  from  the  train  set.  Indeed,  we  can  add  some  new  random  points 
to  the  train  set  to  make  the  size  of  the  train  set  of  constant  cardinality. 

Notice  that  for  the  unchanged  part  of  the  train  set,  we  can  still  apply  the  saturation 
assumption  based  algorithm  to  save  working  load,  and  thus  combine  the  two  ideas. 

The  algorithm  3  is  the  pseudo-code  of  the  adaptively  enriching  greedy  algorithm. 

Remark  3.5.  The  algorithm  can  further  be  modified  in  the  way  that  any  new  parameter  point 
in  retrain  is  subject  to  some  optimization  procedure.  For  instance,  one  can  apply  a  random  walk 
with  decreasing  step  size  and  only  accept  a  new  state  if  the  error  estimator  is  increased.  Or, 
a  more  complex  procedure  such  as  the  simulated  annealing  algorithm  can  be  applied.  Such  a 
procedure  will  additionally  (at  some  cost  though)  increase  the  quality  of  each  added  parameter 
point. 

4  Application  to  the  empirical  interpolation  method 

In  the  following  we  study  how  the  previous  ideas  can  be  used  to  improve  the  greedy  algorithm 
incorporated  in  the  empirical  interpolation  method  (EIM). 


Input:  M:  the  number  of  sample  points  in  each  round  of  searching, 

Nsc:  the  number  of  sample  points  to  pass  the  safety  check, 

Csa,  and  a  tolerance  tol. 

Output:  Sn  and  Wjy 
1:  Set  iVsafe  =  ceil(iVsc/M); 

2:  Generate  an  initial  train  set  Strain  with  M  parameter  samples  (randomly,  or  do  your 
best); 

3:  Choose  an  initial  parameter  value  /i1  £  Strain  and  set  S\  =  {/r1}  and  N  =  1; 

4:  Set  a  vector  i]saved  with  f]saved{n)  =  oo  for  all  p  ^  Strain] 

5:  Compute  set  W\  =  {^(/r1)},  set  safe  =  0,  errortmpmax  =  2  *  tol ; 

6:  while  ( errortmpmax  >  tol  or  safe  <  iVsafe)  do 
7:  Off  OT  tmpmax  =  0; 

8:  for  all  p  £  H  do 

9:  if  Csarisaved{.H)  >  OffOftmpmax  then 

10:  Compute  7/(/u;  WN )  ,  and  let  r]saved( v)  =  ??(/h  Wjv); 

11:  if  Vsaved(^)  >  errortmpmax  then 

12:  OJ'  fOr imprnrix  —  V saved (/■*')  •  and  let  [Amax  —  /h 

13:  end  if 

14:  if  risavedin)  <  tol  then 

15:  flag  /r;  //  all  flagged  parameters  will  be  removed  later 

16:  end  if 

17:  end  if 

18:  end  for 

19:  if  errortmpmax  >  tol  then 

20:  Choose  [an+1  =  Umax,  set  iSjv+i  =  Sn  U  /an+1] 

21:  Compute  r(/rAr+1),  set  Wn+ i  =  Wn  U  {v(fiN+1)}- 

22:  Discard  all  flagged  parameters  from  Strain  and  their  corresponding  saved  error 

estimation  in  iqsaved ; 

23:  Generate  M  —  sizeof(Htrajn)  new  samples,  add  them  into  Strain  such  that 

sizeof (Strain)  =  M;  set  rjsaved  of  all  new  points  to  oo; 

24:  N  4—  N  +  1; 

25:  safe  =0; 

26:  else 

27:  safe  =  safe  +1; 

28:  Discard  Strain ,  generate  M  new  parameters  to  form  Strain  and  set  ft  saved  to  be  oc  for 

all  (A  £  Strain] 

29:  end  if 

30:  end  while 

Algorithm  3:  An  Adaptively  Enriching  Greedy  Algorithm 


9 


4.1  Saturation  assumption 

We  first  test  the  saturation  assumption  based  algorithm  for  two  test  problems  with  low  dimen¬ 
sional  parameter  spaces. 


Test  4.1.  Consider  the  complex-valued  function 

Akx  _  -1 

Fi(x-k)  = - 

x 

where  x  G  12  =  (0,  2]  and  k  G  V  =  [1, 25].  The  interval  12  is  divided  into  a  equidistant  point  set 
of  cardinality  15’000  points  to  build  12/j  where  the  error  ||.Fi(x;/li)  —  ^)\\L°°(nh)  is 

computed.  For  the  standard  greedy  algorithm,  the  train  set  retrain  consists  of  5’000  equidistant 
points  in  V. 


Test  4.2.  As  a  second  and  slightly  more  complicated  example,  consider  the  complex-valued 
function 

^2(x;p)=e^x 

where  the  directional  unit  vector  k  is  given  by 

k  =  —  (sin#  cos  ip.  sin  #  sin  <^,  cos#)2" 


and  fi  =  ( k ,  0,ip)  G  V  =  [1, 5]  x  [0, 7r]  x  [0,  27r] .  As  domain  12  we  take  a  unit  sphere.  For  practical 
purpose,  we  take  a  polyhedral  approximation  to  the  sphere,  as  illustrated  in  Figure  1,  and  the 
discrete  number  of  points  12 h  (where  again  the  error  ||JF2(x;  p,)  —  Xat(J:2)(x;  p)|| L°°(nh)  is  com¬ 
puted)  consists  of  the  three  Gauss  points  on  each  triangle.  For  the  standard  greedy  algorithm, 
the  train  set  Strain  consists  of  a  rectangular  grid  of  50  x  50  x  50  points.  In  computational 


Figure  1:  Discrete  surface  for  the  unit  sphere. 

electromagnetics,  this  function  is  widely  used  since  p.F2(x;/i)  represents  a  plane  wave  with 
wave  direction  k  and  polarization  p  _L  k  impinging  onto  the  sphere.  See  [8] . 

In  Figure  2  we  show  the  evolution  of  the  average  and  maximum  value  of 


C(N) 


#(/A  WN) 

r](p,WN-  i) 


(4.15) 


over  retrain  along  with  the  standard  greedy  sampling  process  Algorithm  1.  This  quantity  is  an 
indication  of  Csa.  We  observe  that  assuming,  in  this  particular  case,  that  Csa  =  2  is  a  safe 
choice,  for  both  cases.  For  this  particular  choice  of  Csa ,  we  illustrate  in  Figure  3  the  savings 
in  the  loop  over  Etrain  at  each  iteration  of  the  standard  greedy  sampling.  Indeed,  the  result 
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Figure  2:  Evolution  of  the  quantity  (4.15)  along  the  greedy  sampling  for  T\  (left)  and  J~2 
(right). 


indicates,  that  at  each  step  N,  the  percentage  of  the  work  that  needs  to  be  done  by  using  the 
saturation  assumption  compared  to  using  the  standard  greedy  algorithm  and  thus  compares 
the  different  workloads,  at  each  loop  over  Strain  of  Algorithm  1  and  2.  One  observes  that, 
while  for  the  first  example  the  improvement  is  already  significant,  one  achieves  an  average 
(over  the  different  values  of  N )  saving  of  workload  of  about  50%  (dotted  red  line)  for  the 
second  example. 


Figure  3:  Percentage  of  work  at  each  step  N,  using  the  Saturation  Assumption  based  greedy 
algorithm,  compared  to  the  workload  using  the  standard  greedy  algorithm,  for  T\  (left)  and 
(right). 


4.2  Adaptively  enriching  greedy  algorithm 

Let  us  also  test  the  adaptively  enriching  greedy  algorithm  (for  convenience  denoted  by  AEGA) 
first  for  the  previously  introduced  function  J-2,  and  then  for  two  test  problems  with  high 
dimensional  parameter  spaces. 

Considering  J-2.  we  set  M  =  l’OOO,  5’000,  25’000,  Nsc  =  125’000  and  Csc  =  2.  In  Figure  4 
we  plot  the  convergence  of  the  new  algorithm  (in  red  solid  lines)  and  the  corresponding  error 
over  the  large  control  set  Strain  (hr  dotted  lines)  consisting  of  125’000  equidistant  points.  For 
comparison,  we  illustrate  the  convergence  of  the  standard  greedy  algorithm  using  the  train  set 
Strain  (in  dashed  lines). 

We  observe  that  as  we  increase  the  value  of  M,  the  convergence  error  provided  by  the 
AEGA  and  the  error  of  the  AEGA  over  the  larger  control  set  become  (not  surprisingly)  closer 
and  closer  and  converge  to  the  error  provided  by  the  standard  EIM  using  training  set  retrain- 
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N  N  N 

Figure  4:  Convergence  behavior  of  the  adaptively  enriching  greedy  algorithm  for  J- 2  with 
M  =  l’OOO  (left),  M  =  5’000  (middle)  and  M  =  25’000  (right). 


Figure  5:  Evolution  of  the  number  of  points  where  the  accuracy  is 
checked  of  the  adaptively  enriching  greedy  algorithm  for  different 
values  of  M  (Safety  check  not  included). 


Figure  5  shows  the  evolution  of  the  number  of  points  which  were  part  of  the  train  set  during 
new  greedy  algorithm  for  all  values  of  M.  We  observe  that  in  all  three  cases  the  error  was 
also  tested  on  at  least  125 ’000  different  points  over  the  iterations  of  the  algorithm,  but  of  low 
number  M  at  each  iteration. 

In  Fig.  6,  we  present,  for  all  values  of  M,  two  quantities.  The  first  one  consists  of  the 
percentage  of  work  (w.r.t.  M),  at  each  step  N ,  that  needs  to  be  effected  and  cannot  be 
skipped  by  using  the  saturation  assumption.  The  second  one  consists  of  the  percentage  of 
points  (w.r.t.  M )  that  remain  in  the  train  set  after  each  step  N.  One  can  observe  that  at  the 
end,  almost  all  points  in  the  train  set  are  withdrawn  (and  thus  the  corresponding  errors  need 
to  be  computed).  During  a  long  time,  we  observe  that  the  algorithm  works  with  the  initial 
train  set  until  the  error  tolerance  is  satisfied  on  those  points  before  they  are  withdrawn.  In 
consequence,  the  accuracy  need  only  to  be  checked  for  a  low  percentage  of  points  of  the  trial 
set  using  the  saturation  assumption.  Towards  the  end,  the  number  of  points  that  remain  in 
the  train  set  after  each  iteration  decreases  and  consequently  for  each  new  sample  point  in  the 
train  set  the  saturation  assumption  cannot  be  used. 

Remark  4.1.  Algorithm  3  is  subject  to  randomness.  However,  we  plot  only  one  realization  of 
the  adaptively  enriching  greedy  algorithm.  Due  to  the  presence  of  a  lot  of  repeated  randomness 
(each  newly  generated  parameter  value  o/H  is  random)  these  illustration  are  still  meaningful. 


Test  4.3.  Introduce  the  following  real-valued  function 

-F3(x;  fj,)  =  sin(27r/ii(xi  -  /i2))  sin(27r/i3(x2  -  /m))  sin(47r/r5(xi  -  /a6))  sin(47r^7(x2  -  ia8)) 
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Figure  6:  Percentage  of  work  (effected  at  each  step  N)  w.r.t.  to  the  total  number  of  points 
M  and  of  the  number  of  points  remained  in  the  train  set  (at  each  step  N)  of  the  adaptively 
enriching  greedy  algorithm  for  J- 2  and  different  values  of  M  =  l’OOO,  5’000, 25’000. 


with  x  =  (xi,x2)  eSl  =  [0,  l]2  and  m  i,M3  G  [0,2],  M2,  M4, M6,  Ms  G  [0,1],  M5,M7  G  [1,2].  The 
domain  11  is  divided  into  a  grid  of  100  x  100  equidistant  points  to  build  11  /, .  Recall  that  in  the 
implementation,  is  used  to  compute  the  norm  ||  •  ||Loom). 

In  Fig.  7  and  8  we  plot  the  convergence  behaviour  of  the  adaptive  enriching  greedy  algo¬ 
rithm  for  the  function  T3  and  tol  =  10~4.  The  value  of  Nsc  is  set  for  all  different  choices  of 
M  =  100,  M  =  l’OOO  and  M  =  lO’OOO  equal  to  Nsc  =  lOO’OOO.  Figure  8  is  a  zoom  of  Figure  7 
towards  the  end  of  the  convergence  to  highlight  how  the  safety  check  acts.  We  observe  that  in 
the  case  of  M  =  lO’OOO  the  safety  check  is  passed  in  only  a  few  attempts  whereas  for  M  =  100 
the  safety  check  fails  34  times  until  finally  successful.  This  means  that  during  each  safety  check 
there  was  at  least  one  parameter  value  where  the  interpolation  error  was  above  the  tolerance. 
Finally,  passing  the  safety  check  means  that  the  interpolation  error  on  lOO’OOO  random  sample 
points  was  below  the  prescribed  tolerance,  in  all  three  cases. 

In  Fig.  9,  the  accumulated  number  of  points  belonging  to  the  train  set  H  is  illustrated.  We 
observe  an  increase  of  this  quantity  towards  the  end  where  the  safety  check  is  active.  During 
this  period  all  parameter  points  are  withdrawn  at  each  iteration,  leading  to  the  increase. 


Test  4.4.  Finally  we  consider  the  following  real-valued  function 


JF^x;  (j.)  =  1  +  exp  — 


(aq-Mi)2  (x2-^2)2 


1  +  exp  - 


M9  M10 

{x\  -  M s)2  O2  -  Me)2 


1  +  exp  - 


(xi  -  M3)2  (x2  -  M4)2 


Ml3 


M14 


1  +  exp  - 


Mil  M 12 

(xi  ~  M7)2  (x2  -  Ms)2 


Mis 


Mi6 


with  x  =  (xi,X2)  G  D  =  [0, 1] 2  and  mi,---,Ms  G  [0.3, 0.7],  M9,---,Mi6  G  [0.01,0.05].  The 
domain  D  =  [0,  l]2  is  divided  into  a  grid  of  100  x  100  equidistant  points  to  build  Qh. 

In  Fig.  10,  the  convergence  behavior  of  the  adaptive  enriching  greedy  algorithm  for  the 
function  J- 4  and  tol  =  10~4  is  plotted.  The  value  of  Nsc  is  set  for  all  different  choices  of 
M  =  lO’OOO,  M  =  lOO’OOO  and  M  =  l’OOO’OOO  equal  to  Nsc  =  lO’OOO’OOO.  Fig.  11  is 
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Figure  7:  Convergence  behavior  of  the 
adaptive  enriching  greedy  algorithm  for  J-\ 3. 


- M  =  1 00 


Figure  8:  Convergence  behavior  of  the 
adaptive  enriching  greedy  algorithm  for  J- 3 
zoomed  in  towards  the  end  where  the  safety 
check  is  active. 


Figure  9:  Evolution  of  the  number  of  points  where  the  accuracy  is 
checked  of  the  adaptively  enriching  greedy  algorithm  for  different 
values  of  M  (Safety  check  included). 


again  a  zoom  of  Fig.  10  towards  the  end  of  the  convergence.  We  observe  a  similar  behavior 
as  in  the  previous  case.  Note  that  in  all  three  cases,  the  safety  check  is  passed  and  the 
tolerance  is  satisfied  for  lO’OOO’OOO  subsequent  parameter  points.  A  bit  surprisingly,  both 
cases  of  M  =  lO’OOO  and  M  =  lOO’OOO  results  in  the  same  number  of  added  modes  N. 

5  Application  to  the  Reduced  Basis  Method 

In  this  section  we  apply  the  improved  greedy  algorithms  in  the  context  of  the  Reduced  Basis 
Method  (RBM).  As  in  the  last  section,  we  first  test  the  benefit  of  the  saturation  assumption 
for  some  low  dimensional  parametric  problems,  and  then  proceed  to  test  a  problem  with  15 
parameters  using  the  adaptively  enriching  greedy  algorithm. 

5.1  Saturation  Assumption 

Before  performing  the  numerical  test,  we  shall  show  that  for  some  variational  problems,  the 
saturation  assumption  is  satisfied  with  Csa  =  1  if  the  error  is  measured  in  the  intrinsic  energy 
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Figure  10:  Convergence  behavior  of  the 
adaptive  enriching  greedy  algorithm  for  J-\. 


Figure  11:  Convergence  behavior  of  the 
adaptive  enriching  greedy  algorithm  for  J- 4 
zoomed  in  towards  the  end  where  the  safety 
check  is  active. 


norm. 

Suppose  the  variational  problem  is  based  on  a  minimization  principle, 

u  =  argmin^gjf  J(v),  (5.16) 

where  J(v )  is  an  energy  functional.  Then  the  finite  element  solution  u*e  on  X?e  C  X  is 

ufe  =  avgmmveXfe  J(v).  (5-17) 

Similarly,  the  reduced  basis  solution  urJ}  on  X $  C  X^e  C  X  is 

urN  =  ar  gimnveXrbJ(v).  (5.18) 

The  error  between  u*e  and  ux  can  be  measured  by  the  nonnegative  quantity 

J{ur£)-  (5.19) 

Since  Wn  C  Wn+i  and  thus  X x  C  Xrx+l,  we  have  J(urx+1)  <  J(u\ ^)  and  in  consequence 
there  holds 

j(t$+i)  -  J(ufe)  <  J(urj$)  -  J(uf£).  (5.20) 

Observe  that  if  the  error  is  measured  exactly  as  J(u\ y)  —  J(t/e),  the  saturation  assumption  is 
satisfied  with  Csa  =  1. 

Let  us  give  an  example  of  the  above  minimization  principle  based  problem.  Consider  the 
symmetric  coercive  elliptic  problem, 

(  — V  •  (a(/i)Vri)  =  /  in  Q, 

<  u  =  0  on  r^i,  (5.21) 

[  a(g)V«'  n  =  g  on  Tjv, 

where  T £>  and  T jy  are  the  Dirichlet  and  Neumann  boundaries  of  XI,  and  T d  U  T jy  =  <9fL 

For  simplicity,  we  assume  7^  0.  The  functions  /  and  g  are  L2-functions  on  12  and  F jy 
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respectively.  The  parameter  dependent  diffusion  coefficients  a(/i)  are  always  positive.  Let 
X  =  H  )-,(&,)  :=  {«£  :  u|rD  =  0}.  Then  its  variational  formulation  is 

a(u,v;  /i)  =  f(v),  Vv  £  X,  (5.22) 


where 

a(u ,  v ;  /^) 

The  energy  functional  is 


a(fi)XuXvdx, 


f(v ) 


/  fvdx + 

Jn 


gvds. 


J(v )  =  ^ll(a(/^))1/2Vt>||o,n  -  /(«)* 


and  we  have 

J«)  -  =  1  (||a1/2V«S||^  -  ||aVW'||iS,n)  -  /(u#  -  «/') 

=  i||c,V2V(«^  -  a/e)||S,n  +  J  <*V«/e  •  V(<  -  ,/e)<ix  -  /(»$  -  t/') 
=  ^||ci1/2V(ttS-a/e)l|§,a. 

In  the  last  step,  Galerkin  orthogonality 


f  aSJu^e  ■  Xiudx  =  f(w )  V 

Jci 


w 


e  x/e 


is  used.  The  above  analysis  is  standard,  see  [13]. 

If  we  measure  the  error  by  this  intrinsic  energy  norm  ||a(/i)1/2Vu||o,f2,  the  error  satisfies 
the  saturation  assumption  with  Csa  =  1, 


a(/i)1/2V {urM  -  u/e)||0,n  <  || a(/i)1/2V(u^  -  u/e)||0,n  for  0  <  TV  <  M.  (5.23) 


Unfortunately,  even  for  the  above  model  problem,  the  a  posteriori  error  estimator  used  in  the 
reduced  basis  method  is  not  equivalent  to  the  energy  norm  of  the  error  exactly. 

For  the  problem  (5.21),  we  choose  the  underlying  norm  of  X  and  X^e  to  be  the  H1- semi- 

norm  ||r>||x  =  \JJq (Vu)2^.  Notice,  due  to  the  fact  that  the  dual  norm  of  the  X-U-norni 

needs  to  be  computed,  involving  an  inverse  of  the  matrix  associated  with  the  X4e-norm,  this 
Xfe  -norm  cannot  be  parameter  dependent.  Otherwise,  we  must  invert  a  matrix  of  a  size 
comparable  to  the  finite  element  space  every  time  for  a  new  parameter  in  the  computation  of 
error  estimator.  This  would  clearly  result  in  an  unacceptable  cost. 

The  functional  based  error  estimator  for  (5.21)  is  defined  as 


wn) 


II/II.y'IK-;m)IL-' 
M\fiuTN^))\  ’ 


(5.24) 


For  this  simple  problem  and  this  specific  choice  of  norm,  the  stability  constant  is  /3(/x)  = 
minxen{a(AA)}-  Note  that 

f{v)  =  a(v,  v;  fi)  =  ||a(/x)1/2Vu||^n,  Vu  G  X. 


The  error  estimator  rj(n;  Wn)  we  used  here  is  in  fact  an  estimation  of  the  relative  error  mea¬ 
sured  by  the  square  of  the  intrinsic  energy  norm.  In  principle,  since  we  already  proved  the 
error  measured  in  energy  norm  should  be  monotonically  decreasing  (or  more  precisely,  non¬ 
increasing),  the  constant  Csa  can  be  chosen  to  be  1.  However,  if  we  examine  the  error  estimator 
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r/(/x;  Wat)  defined  in  (5.24)  carefully,  we  find  that  for  a  fixed  parameter  /x,  the  values  of  ||/||.y' 
and  /3(/x)  are  fixed,  the  change  of  the  value  of  \f(ur^(fi))\  is  quite  small  if  the  approximation 
ur^(fi))  is  already  in  the  asymptotically  region.  The  only  problematic  term  is  || ?’(•;  /i)||x'-  This 
term  is  measured  in  a  dual  norm  of  a  parameter  independent  norm  (the  H 1-semi-norm) ,  not 
in  the  dual  norm  of  the  intrinsic  energy  norm  ||a(/^)1/2Vu||o,n.  It  is  easy  to  see  that 


a1/2Ve 


||  2  _ 

Il0,fi  — 


/(e)  < 


x'||e||x  < 


X' 


/?(/*) 


Thus,  the  error  estimator  r/(/x;  Wjsr)  is  only  an  upper  bound  of  the  the  relative  error  measured 
by  the  square  of  the  intrinsic  energy  norm.  When  the  variation  of  a  with  respect  to  /x  is 
large,  the  difference  between  the  H 1  -semi-norm  and  the  energy  norm  ||a1/,2Vu||o,n  may  be 
large  and  we  may  find  the  error  estimator  is  not  monotonically  decreasing  as  the  number  of 
basis  functions  increasing. 

In  Test  5.1  below,  we  use  a  moderate  variation  of  a  £  [1/10,10]  and  we  observe  that 
the  saturation  assumption  is  satisfied  with  Csa  =  1.  In  Tests  5.2  and  5.3,  a  wider  range  of 
a  G  [1/00,100]  is  used.  For  the  corresponding  error  estimator,  the  saturation  assumption  is 
not  satisfied  with  Csa  =  1  but  for  a  larger  Csa. 


Test  5.1.  In  this  example,  we  will  show  that  for  a  simple  coercive  elliptic  problem,  the 
saturation  assumption  is  satisfied  numerically  with  Csa  =  1  for  some  type  of  error  estimators. 

Consider  the  following  thermal  block  problem,  which  is  a  special  case  of  (5.21),  see  also 
[14].  Let  17  =  (0,  l)2,  and 


'  — V  •  (aVu) 

=  1 

in 

n, 

u 

=  0 

on 

Ttop  =  {iG  (0, 1  ),y  =  1}, 

aVu  ■  n 

=  0 

on 

F side  =  n  =  0  and  x  =  l,y  G  (0,1)} 

aVu  ■  n 

=  1 

on 

Fbase  =  n  G  (0, 1  ),y  =  0}, 

(5.25) 


where  a  =  lO2^”1  in  R\  =  (0,0. 5)2  and  a  =  1  in  Rrest  =  17\(0,0.5)2.  We  choose  the  one- 
dimensional  parameter  domain  V  of  /x  to  be  [0, 1],  which  corresponds  to  a  £  [1/10, 10]  in  R.\ . 
Note  that  in  this  particular  case  the  vector  of  parameters  /x  is  indeed  a  scaler  parameter  fi. 
The  variational  problem  is  given  in  (5.22).  The  output  functional  is  chosen  to  be  s(u )  =  f{u). 

Let  T  be  a  uniform  mesh  on  17  with  80’401  of  nodes  (degrees  of  freedom),  and  P\(K)  be 
the  space  of  linear  polynomials  on  an  element  K  G  T.  We  then  define  our  finite  approximation 
space 

I/e  =  {»Gl:  v\K  G  Pi(K),  \/K  G  T}. 

For  a  given  /x,  the  finite  element  problem  is  seeking  we{^)  G  X*e,  such  that 


a(ufe(ii),v;  n)  =  f(v)  vGl/e. 


(5.26) 


With  a  set  of  N  reduced  basis  elements  Ifby  and  the  corresponding  X$,  and  for  a  given 
parameter  /x,  we  solve  the  following  reduced  basis  approximation  problem:  Seeking  G 

X(/,  such  that 

a(ur^(/j,),v,n)  =  f(v)  VrGljJ.  (5.27) 

We  choose  101  equidistance  sample  points  in  V  =  [0, 1],  i.e.,  Strain  =  {^,7  =  0,1, 2, •••  ,100}, 
a  standard  reduced  basis  greedy  algorithm  with  the  error  estimator  defined  in  (5.24)  being  used. 
The  first  parameter  /x1  is  chosen  to  be  0.5,  that  is,  S\  =  {0.5}.  The  greedy  algorithm  chooses 
the  2nd,  3rd,  4th,  and  5th  parameters  as  {0,1.0,0.15,0.81},  so  S$  =  {0.5,0,1,0.15,0.81}.  We 
compute  the  reduced  basis  set  Wn  and  the  reduced  basis  spaces  X$  corresponding  to  Sn, 
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N  =  1,  ■  ■  ■  ,  5.  Then,  for  all  points  /i  in  Etrain,  Wn),  N  =  1,  •  •  •  ,  5  is  computed.  Figure  12 
shows  the  plots  of  rj  for  each  points  in  Strain  with  number  of  reduced  basis  =  1,  ■  ■  •  ,5.  Clearly, 
we  see  that  for  each  point  G  Strain,  viv,  Ws)  <  r/(/x;  W4)  <  r/(/x;  W3)  <  r?(/x;  W2)  <  77(7x5  Wi). 
For  this  problem,  the  saturation  assumption  is  clearly  satisfied  for  Csa  =  1  in  the  first  several 
steps. 


Figure  12:  A  verification  of  Saturation  Assumption  for  a  symmetric  positive  definite  problem 
with  1  parameter. 


Test  5.2.  We  now  test  the  potential  for  cost  savings  when  the  saturation  assumption  based 
algorithm  is  used  in  connection  with  the  reduced  basis  method. 

For  equation  (5.25),  we  decompose  the  domain  12  into  3  subdomains:  R\  =  (0, 0.5)  x  (0.5, 1), 
R2  =  (0.5, 1)  x  (0,  0.5),  and  R3  =  (0,  l)2\(i?i  U  R2).  The  diffusion  constant  a  is  set  to  be 


(  a,  =  1002«_1,  xeRi,i  =  1,2, 
1  «3  =  1,  X  €  R3, 


where  /x  =  (/r  1,7x2)  £  [0,  l]2.  The  domain  of  a^i  =  1,  2  is  set  to  [1/100, 100].  The  bilinear  form 
then  becomes 

a(u,  v]  fi)  =  1002Mi_1  f  Vu-'Vvdx+  f  Vu-Vvdx.  (5.28) 

JRi  Jr3 


All  other  forms  and  spaces  are  identical  to  the  ones  of  Test  5.1.  Further,  let  Woo  =  {0, 1, 2,  •  •  •  ,  100}, 
the  train  set  is  given  by 

Strain  =  {( x{i),y{j ))  :  x{i)  =  =  1^0 ,  for  i  G  Woo,  j  e  Woo}- 

We  set  the  tolerance  to  be  10-3  and  use  the  error  estimator  defined  in  (5.24). 

Both  the  standard  greedy  algorithm  and  the  saturation  assumption  based  greedy  algorithm 
requires  20  reduced  bases  to  reduce  the  estimated  error  to  less  than  the  tolerance.  For  this 
problem,  if  the  error  is  measured  in  the  intrinsic  energy  norm,  we  can  choose  Csa  =  1  as 
indicated  above.  Due  to  the  inaccuracy  of  the  error  estimator,  the  saturation  assumption 
based  algorithm  chooses  a  slightly  different  set  of  Sjv-  If  we  choose  Csa  =  1.1  slightly  larger, 
we  get  however  the  same  set  Sjy  as  the  standard  greedy  algorithm. 
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See  Fig.  13  for  a  comparisons  of  the  workloads  using  the  standard  algorithm  and  the 
new  approach  based  on  the  saturation  assumption  with  Csa  =  1  and  Csa  =  1.1,  respectively. 
The  mean  percentage  is  computed  as  X^i(PercentaSe  at  step  *)/AL  The  mean  percentages 
of  Csa  =  1  and  Csa  =  1.1  is  about  32%  and  34%  respectively.  Since  the  computational  cost  for 
each  evaluation  of  r/(/x;  N)  is  of  0(N3),  the  average  percentages  do  not  represent  the  average 
time  saving,  and  only  give  a  sense  of  the  saving  of  the  workloads  at  each  step. 

In  Fig.  14,  we  present  the  selections  of  the  sample  points  Sjy  by  the  standard  algorithm 
and  the  Saturation  Assumption  based  greedy  algorithm  with  Csa  =  1.  Many  sample  points 
are  identical.  In  the  case  Csa  =  1.1,  the  samples  Sjy  are  identical  to  the  ones  of  the  standard 
algorithm. 


Figure  13:  Percentage  of  work  at  each  step  N  Figure  14:  Sjy  obtained  by  standard  and  sat- 
using  saturation  assumption  based  greedy  al-  uration  assumption  based  greedy  algorithms 
gorithm  with  Csa  =  1,  compared  to  the  work-  for  Test  5.2. 
load  using  the  standard  greedy  algorithm  for 
Test  5.2. 

Test  5.3.  We  continue  and  test  a  problem  with  3  parameters.  For  (5.25),  we  decompose  the 
domain  fl  into  4  subdomains:  Rk  =  (%p,  |)  x  (^-,  |),  for  i  =  1, 2,  j  =  1, 2,  and  k  =  2(i  —  l)+j. 
The  diffusion  constant  a  is  set  to  be 

_  f  ak  =  lOO2^-1  x  G  Rk,k  =  1,2,3, 

}  CC4  =  1  X  E  i?4, 

where  n  =  (fii,  /j,2,  fa)  G  [0,  l]3.  The  bilinear  form  is 

1002/ifc_1  /  Vu  -Vvdx  +  /  Vu  •  Vvdx,  (5.29) 

J  Rk  ^  R& 

All  other  forms  and  spaces  are  identical  to  the  ones  of  Test  5.1.  We  again  choose  the  output 
functional  based  error  estimator  as  Test  4.1  and  the  tolerance  is  set  to  be  10~3.  Let  IV50  = 
{0, 1, 2,  •  •  •  ,  50},  the  train  set  is  given  by 

Strain  =  {(a;(*)  >  y{j)  >  z{^))  ■  x{l)  =  gQ  ,y{j)  =  ,  z{k)  =  G  A50,  j  G  -/V50,  k  G  AI50}. 

The  standard  greedy  algorithm  needs  24  reduced  bases  to  reduce  the  estimated  error  less  than 
the  tolerance.  If  Csa  is  chosen  to  be  1,  25  reduced  bases  are  needed  to  reduce  the  estimated  error 
less  than  the  tolerance.  The  set  Sn  obtained  by  the  Saturation  Assumption  based  algorithm 


3 

a{u,v,n)  = 

k=  l 
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with  Csa  =  1  is  also  different  from  the  standard  algorithm.  As  discussed  before,  this  is  mainly 
caused  by  the  inaccuracy  of  the  error  estimator.  If  we  choose  Csa  =  3,  we  obtain  the  same 
sample  points  Sjv  as  the  standard  greedy  algorithm.  See  Figure  15  for  the  comparisons  of  the 
workloads  using  the  standard  algorithm  and  the  saturation  assumption  based  algorithm  with 
Csa  =  1  and  Csa  =  3,  respectively.  The  mean  percentages  of  workload  for  Csa  =  1  and  Csa  =  3 
are  21.6%  and  33.7%,  respectively. 


N 


Figure  15:  Percentage  of  work  at  each  step  N 
using  saturation  assumption  based  greedy  al¬ 
gorithm  with  Csa  =  1  and  Csa  =  3,  compared 
to  the  work  load  using  the  standard  greedy 
algorithm  for  Test  5.3. 


Figure  16:  Convergence  behavior  of  the 
adaptively  enriching  greedy  algorithm  for 
Test  5.4  with  M  =  100,  500,  and  l’OOO. 


Remark  5.1.  For  the  type  of  compliance  problem  discussed  in  Tests  5.1,  5.2  and  5.3,  other 
types  of  error  estimator  are  suggested  in  [14]: 


rf(» ,  WN) 


H-;/*)ll  (Xfey 

/?/*(/*)V2|<(/*)  I 


and  3S{^WN) 


IK-;/*)ll  2(Xf°y 

Pfe(v)\ Sjv(m)I  ' 


(5.30) 


As  discussed  above,  the  most  important  term  in  the  error  estimator  of  the  Saturation  Assump¬ 
tion  is  the  dual  norm  of  the  residual  n)\\txfe)' ■  F°r  the  error  estimator  r]e(n]Wisr),  the 
behavior  is  similar  to  that  of  ri(pi]  Wn).  For  the  error  estimator  r/s(fi]  Wjy),  the  dual  norm  of 
the  residual  is  squared.  The  dual  norm  is  computed  with  respect  to  a  parameter  independent 
reference  norm.  The  square  makes  the  difference  between  the  dual  norm  based  on  the  intrinsic 
energy  norm  and  on  the  reference  norm  larger.  Normally,  if  the  error  estimator  Wn)  is 
used,  we  need  a  more  conservative  Csa.  Numerical  tests  show  that  even  if  Csa  =  20  is  set  for 
Test  5.3,  the  workload  of  the  saturation  assumption  based  algorithm  is  still  only  about  45%  (on 
average)  of  the  standard  greedy  algorithm. 


5.2  Adaptively  enriching  greedy  algorithm 


Test  5.4  We  test  the  adaptively  enriching  greedy  algorithm  for  the  reduced  basis  method  for 
a  problem  with  15  parameters. 

For  (5.25),  we  decompose  the  domain  12  into  16  subdomains:  Rk  =  (^,  f)  x  f),  for 
i  =  1,2,  3, 4,  j  =  1,  2,  3, 4,  and  k  =  4 (i  —  1)  +  j.  The  diffusion  constant  a  is  set  to  be 


x  €  Rk, 
x  E  R\q. 


k  =  1,2,  -  -  •  ,15, 


a  = 


046 


52/ifc-l 

1, 
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where  /x  =  (/xi ,  ^2,  ■  ■  •  ,  /X15)  E  [0,  l]15.  The  domain  of  ccfe,  fc  =  1,2,...,  15,  is  given  by  [1/5,  5]. 
The  bilinear  form  then  consists  of 


15 

a(u,v;n )  =  ^52m 
fc=i 


fe-i 


\7u  ■  Vvdx  + 


\7u  ■  Vvdx. 


'  Rk 


'  Rl6 


(5.31) 


All  other  forms  and  spaces  are  identical  to  the  ones  of  Test  5.1.  Due  to  the  many  jumps  of  the 
coefficients  along  the  interfaces  of  the  subdomains,  the  solution  space  of  this  problem  is  very 
rich.  We  set  Csa  =  1,  tol  =  0.05,  Nsc  =  lO’OOO.  Since  there  is  a  “safety  check”  step  to  ensure 
the  quality  of  the  reduced  bases,  we  should  not  worry  that  the  choice  of  constant  Csa  is  too 
aggressive.  We  test  three  cases:  M  =  100,  M  =  500,  and  M  =  l’OOO.  The  convergence  for 
one  realization  is  plotted  in  Fig.  16.  The  number  of  reduced  basis  for  M  =  100  is  52,  and  for 
the  other  two  cases  is  50.  This  suggests  us  that  a  bigger  M  will  lead  to  a  smaller  number  of 
the  bases.  The  percentage  of  work  (effected  at  each  step  N)  with  respect  to  the  total  number 
of  points  M  and  of  the  number  of  points  remained  in  the  train  set  (at  each  step  N)  of  the 
adaptively  enriching  greedy  algorithm  for  different  values  of  M  =  100,  500,  and  1000  are  shown 
in  17.  At  the  beginning  stage,  the  estimated  errors  are  larger  than  tolerance  for  almost  all 
points  in  the  train  set,  when  the  RB  space  is  rich  enough,  more  and  more  points  are  removed, 
and  eventually,  almost  all  points  in  the  trained  set  are  removed  in  later  stages.  For  the  number 
of  the  points  where  the  error  estimators  are  computed,  that  is,  the  saturation  assumption  part, 
it  behaves  like  the  Algorithm  2  with  a  fixed  train  set  since  the  train  set  barely  changed  in  the 
beginning.  In  the  later  stage,  since  almost  all  points  are  new  points,  the  percentage  of  the 
number  of  the  points  where  the  error  estimators  are  computed  is  close  to  100%. 


N  (M=500) 


N  (M=1000) 


Figure  17:  Percentage  of  work  (effected  at  each  step  N)  w.r.t.  the  total  number  of  points 
M  and  of  the  number  of  points  remained  in  the  train  set  (at  each  step  N)  of  the  adaptively 
enriching  greedy  algorithm  for  Test  5.4  and  different  values  of  M  =  100,  500,  and  l’OOO. 


Remark  5.2.  For  a  fixed  M  and  provided  the  algorithm  is  performed  several  times,  we  observe 
that  even  though  the  train  set  is  generated  randomly  each  time,  the  numbers  of  the  reduced  bases 
needed  to  reduce  the  estimated  error  to  the  prescribed  tolerance  are  very  similar.  This  means 
that  even  if  we  start  with  a  different  and  very  corse  random  train  set,  the  algorithm  is  quite 
stable  in  the  sense  of  capturing  the  dimensionality  of  the  reduced  basis  space. 
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6  Conclusions 


In  this  paper,  we  propose  two  enhanced  greedy  algorithms  designed  to  improve  sampling  ap¬ 
proaches  for  high  dimensional  parameters  spaces.  We  have  demonstrated  the  efficiency  of  these 
new  techniques  on  the  empirical  interpolation  method  (EIM)  and  the  reduced  basis  method 
(RBM).  Among  other  key  observations,  we  have  documented  the  potential  for  substantial  sav¬ 
ings  over  standard  greedy  algorithm  by  utilization  of  a  simple  saturation  assumption.  Com¬ 
bining  this  with  a  ”  safety  check”  step  guaranteed  adaptively  enriching  greedy  algorithm,  the 
EIM  and  RBM  for  problems  with  a  high  number  of  parameters  are  now  workable  and  more 
robust. 

With  some  possible  modifications,  the  algorithms  developed  here  can  be  applied  to  other 
scenarios  where  a  greedy  algorithm  is  needed,  for  example,  the  successive  constraint  linear 
optimization  method  for  lower  bounds  of  parametric  coercivity  and  inf-sup  constants  [3,  4,  5, 
10]. 
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